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ABSTRACT 

We investigate the gravitational fragmentation of expanding shells driven by HII regions using 
the three-dimensional Lagrangian simulation codes based on the Riemann solver, called Godunov 
smoothed particle hydrodynamics. The ambient gas is assumed to be uniform. In order to attain 
high resolution to resolve the geometrically thin dense shell, we calculate not the whole but a part of 
the shell. We find that perturbations begin to grow earlier than the prediction of the linear analysis 
under the thin-shell approximation. The wavenumber of the most unstable mode is larger than that in 
the thin-shell linear analysis. The development of the gravitational instability is accompanied by the 
significant deformation of the contact discontinuity. These results are consistent with a linear analysis 
presented by Iwasaki et al. that have taken into account the density profile across the thickness and 
approximate shock and contact discontinuity boundary conditions. We derive useful analytic formulae 
for the fragment scale and the epoch when the gravitational instability begins to grow. 
Subject headings: HII regions - hydrodynamics - instabilities - shock waves - stars: formation 



1. INTRODUCTION 

Massive stars strongly disturb the interstellar medium 
(ISM) via emission of ionizing photons, stellar winds, 
and supernova explosions. These processes produce over- 
pressured hot bubbles that expand into ambient inter- 
stellar clouds. At the same time, a shock wave sweeps 
up the ambient clouds into a dense shell. The ex- 
panding shell strongly influences the dynamics of the 
ISM. Especially, it is often supposed to trigger the for- 
mati on of stars through the gr avitational fragmenta- 
tion (lElmegreen fc Lada IflOTTf) . iHosokawa &: Inutsuka I 
(|2005l l2006f ) investigated the dynamical expansion of 
the HII region, the photodissociation region, and the 
swept-up shell, solving the UV and far-UV radiative 
transfer and thermal and chemical processes in the one- 
dimensional (ID) hydrodynamics code. They showed 
that the shell becomes cold and dense enough for the 
gravitational instability (GI) to take place owing to the 
reformation of molecules destructed by far-UV photons. 
Numerous authors have discovered evidences for the star 
forma tion in shell-like mo l ecula, r cloud s around hot bub- 
bles. iChurchwell et al. I (|2006l . I2007D have compiled a 
catalogue of ~ 600 shells from data in Spitze r-GLIMPSE 
survey. Recently, iDeharveng et al.~l ()2010[ ) have inves- 
tigated 102 samples identified as shells on the Spitzer- 
GLIMPSE images at S.O^m. They found that 86% of the 
shells arc associated with ionized gases, or HII regions. 
They obtained statistical properties of the triggered star 
formation, and suggested that more than a quarter of the 
shells may have triggered the formation of massive stars. 
This suggests that the trigger star formation process may 
be important in the massive star formation. 

To understand the triggered star formation, it is im- 
portant to investigate how and when the expanding 
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shell fragments th r ough the GI. lElmegreen I ()1994[ ) and 
iWhitworth et al. 1 ()1994[ ) derived a simple dispersion re- 
lation taking into account the mass accretion and the 
dilution effect of perturbations owing to the expansion. 
They assumed the thin-shell approximation, and neglect 
the bound ary effect of the cont act discontinuity (CD). 
Recently, i Iwasaki et al. I ()2011l hereafter Paper I) in- 
vestigated stability of expanding shells taking into ac- 
count asymmetric density profiles with imposing an ap- 
proximate shock boundary condition on the leading sur- 
face and the CD boundary condition on the trailing sur- 
face. They found that the dispersion relation and eigen- 
function strongly depends on the boundary conditions 
and the degree of asymmetry of the density profile. 

Although many authors have studied fragmentation of 
shells by using linear analysis, it is still uncertain how 
and when shells fragment. This is because the stabil- 
ity analysis of the evolving shell is difficult to perform 
without any approximations. Therefore, time-dependent 
multi-dimensional simulations are crucial to investigate 
it. To resolve propagating thin shells all the time, the 
Eulerian adaptive mesh refinement (AMR) technique in 
the mesh code or the smoothed particle hydrodynamics 
(SPH) have been used. However, since enormous meshes 
or SPH particles are required to resolve the thickness, 
researc hes based on nume rical simulations are less ad- 
vanced. IDale et al. I ()2007f ) investigated the gravitational 
fragmentation of the shell driven by the expansion of the 
HII region by using three-dimensional (3D) SPH taking 
into account the radiative transfer of ionizing photons. 
In their s imulation, i ts thi ckness was not resolved suf- 
ficiently. IDale et all (|2009D have investigated fragmen- 
tation of shells with fixed mass that expand into a hot 
rarefied gas by using two different numerical schemes, the 
Eulerian AMR method (the flash code) and the SPH 
method. They considered shells confined by a constant 
pressure on both surfaces. In the more realistic situa- 
tion where shells are confined by the CD and the shock 
front (SF), multi-dimensional simulations including the 
self-gravity have not been performed yet. 
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In this paper, we perform the 3D simulation to in- 
vestigate the fragmentation process of expanding shells 
through GI. As numerical method, we adopt the 3D SPH 
method. In order to attain high resolution enough to re- 
solve the thickness, we calculate not the whole but a part 
of the shell. The outline of the paper is as follows: in Sec- 
tion [2J we describe our model and simulations in detail. 
Brief review of Paper I is presented in Section [3l The 
results of our simulations are shown in Section |4l In Sec- 
tions [5] and m we present discussions and summaries of 
our paper. 

2. SIMULATIONS 

2.1. Numerical Methods 

We adopt the SPH method (e.g., IMonaghan I I1992D . 
Since SPH is Lagrangian particle method, it is one of 
the best method for problems where a wide low density 
region and geometrically thin shell coexist. Instead of 
additional viscosity to handle shocks, we use the "Go- 
dunov SPH" where the results of the Riemann problem 
are used to calculate the interaction bet ween SPH parti- 
cles (IInutsukall2002| ). The tree method (jl arnes fc Hut I 
119861) is used to calculate self-gravitational force. 

2.1.1. Treatment of Contact Discontinuity 

In the standard SPH method (jMonaghan II1992D . the 
i-th particle has the mass of m^, and its density is given 
by 

Pi ^^mkW{-K^-Xk,hik), (1) 
k 

where W{x, h) is kernel function and we use the Gaussian 
kernel, h is the smoothing length, and hik is an average 
between hi and hk- The equation of motion of the i-th 
particle is given by 
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The standard SPH has a shortcoming wh en calculating 
the C D between a ho t and a cold gas (|Agertz et al. I 
I2007f l. lAgertz et aT~l (|2007l ) reported that the artifi- 
cial tension stabilizes the Kelvin-Helmholtz instability 
with a different d ensities in the standard SPH. In con- 
trast, [Ch^^eF^lT] ([2010) have shown that Godunov SPH 
can describe the Kelvin-Helmholtz instability reasonably 
well. 

In addition. [ Tartakovsky fc Meakin I (|20b5l ) and 
IHu fc Adams I (|2OO60 proposed a new method that treats 
the CD naturally in the context of a multi-phase fluid. 
They modified Equations ([T]) and ^ to treat the CD as 
follows: 

Pi TOi^ W^(xj - Xfe,/ijfe), (3) 
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The term of mi / pi represents the volume of i-th particle 
~ /if. Therefore, if we determine particle mass of hot 



and cold gases so that its smoothing length distributes 
smoothly across the CD, the terms inside the bracket in 
Equation (j4]) is smooth across the CD, leading to much 
better behavior at the CD. Their new method is very 
useful in problems where the position of the CD is known 
in advance as in our modeling. Equation Q satisfies 
the relation of TO^Avi = — rrifcAvfc, suggesting that the 
momentum conservation is guaranteed. Furthermore, we 
apply the Godunov technique to Equation Q as 
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where Pj*^, is the result of the R iemann problem between 
the i-th and the fc-th particles (|Inutsuka Il2002f ). 

2.2. Models 

Massive stars emit ultraviolet photons {hv > 13.6 eV) 
and produce HII regions around them. Here, we con- 
sider a massive star that emits ionizing photons with the 
photon number luminosity Quv [s^^] into the ambient 
gas with the uniform density of pE = totie, where riE 
and m are the number density and the mean mass of the 
ambient gas particle, respectively. 

We construct as simple model as possible to concen- 
trate on the physics of GI of expanding shells. Figure 
[1] illustrates the schematic picture of our model. In this 
paper, we do not calculate the radiative transfer of ion- 
izing photons but introduce hot and cold gases that are 
assumed to evolve keeping their temperatures constant. 
The hot gas is occupied in r < i?cD (the dotted region), 
and the thin shell is in i?cD < r < Rsf, where i?cD and 
RsF are positions of the CD and the SF, respectively. 
The pre-shock ambient cold gas in r > i?sF is assumed 
to be uniform. In our model, the inner boundary is set 
at r = i?b inside the hot bubble (i?b < Rcu)- In the HII 
region, the detailed balance between the recombination 
and the ionization is approximately established all the 
time. Therefore, the interior pressure of the HII region 
is given by 

-Pint = PECii f . (6) 

where i?sT is the Stromgren radius. 
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We assume that the pressure of the HII region is spatially 
constant because of the high temperature. The hot gas 
at r = i?b is pushed by the interior pressure Pint- 

2.3. Calculation Domain 

If we simulate the overall shell, the required number 
of SPH particles, A'^tot, is too large to calculate. There- 
fore, to save -/Vtot, we calculate a part of the shell. The 
schematic picture of the calculation domain is shown in 
Figure m and is designated by \ tan~^{y/x)\ < 9^, and 
I tan~^ (z/x)! < 9t,. Since the solid angle subtended by 
the calculation domain from the center O is given by 

~ A9l for 61b < 1, (8) 
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Fig. 1. — Schematic picture of our model. The hot gas occupying 
the dotted region is embedded by the cold gas. 

the total number of SPH particles can be reduced by a 
factor of ilb/47r ^ S^/t^ compared with the calculation 
of the whole shell. 




Fig. 2. — Schematic picture of the calculation domain. 

2.4. Simulation Units 

In this paper, for convenience, the units of the time, 
length, and mass scales are taken to be 

i?o = (^) = 5.9 Ql^;,,, n,f PC (10) 

and 

Mo = peRI = 5.0 X 10^ Q^/J 49 n-'g/' Mq, (11) 

respectively, where (5uv,49 = Quv/lO"*^ s~^, and nE,3 = 
ng/lO^ cm~3. The dependence of the expansion law on 
the parameters (QuV; it-e) are approximately eliminated 
by using the non-dimensional quantities normalized by 
^0, Ro, and Mq (see Figure 1 in Paper I). 



2.5. Initial and Boundary Conditions 

It is difficult to resolve very thin shell in early evolu- 
tionary phase even if we calculate the part of the shell 
(see Figured]). Thus, we use a grid-based ID hydrody- 
namical code to describe the early evolutionary phase of 
the HII region expansion, and we switch to 3D SPH at 
t = 0.4^0- As a grid-based numerical method, we use the 
ID spherically s ymmetric Lagrangian Godunov method 
(jvan Leer 1119971 ). The innermost mesh at i?b is pushed 
by the interior pressure given by Equation (jS]). When 
the ID simulation reaches t/to = 0.4, the radial profiles 
of physical quantities are used as the initial condition of 
the 3D calculations. We smooth the distribution of the 
sound speed across the CD in the scale of the smoothing 
length. It is assumed that the temperature of the hot gas 
is 64 times as high as that of the cold gas. The specific 
value of the hot gas temperature does not change our 
conclusion as long as it is much larger than the temper- 
ature of the cold gas. 

We calculate the expansions of HII regions around the 
41M0 (the high-mass (HM) model) and 19Mq (the low- 
mass (LM) model) stars that are embedded by the uni- 
form ambient gas of ue = 10^ cm""^. Simulation pa- 
rameters are listed in Table [TJ The temperature of the 
cold gas Tc is assumed to be 10 K. The opening angles 
of simulation domains are set to 9h = 27r/26 and 27r/14 
for the HM and the LM models, respectively, so that the 
calculation domains contain a single wavelength of the 
most unsta ble mode predicted from the thin-shell linear 
analysis by lElmegreen ] (pfll . The mass of SPH par- 
ticles are set so that the initial thickness is resolved by 
five SPH particles. Since the thickness increases with 
time, the relative resolution becomes better as the shell 
expands. In the later phase, the thickness is resolved by 
~ 15 particles. The total number of SPH particles for the 
HM and the LM models are 4.00 x 10^ and 2.26 x 10^, 
respectively. 

Corrugation-type perturbations are put into the shell 
at the initial state. We displace the SPH particles so that 
their densities do not change, or V • Ax^ — 0, where Ax^ 
is the displacement of the z-th particle. The functional 
form of the displacement of the shell is assumed to be 
a — cos{l(p), where </> = ta.n~^ {y / x) . Therefore, the shell 
has the negative displacement at = 0. The SPH parti- 
cles are displaced keeping their velocity, the sound speed, 
and the mass fixed. We concentrate on the evolution of 
a single mode in each calculation. Dependence on I is 
investigated by many runs. 

The moving boundary conditio n at r = Rh(t) is re- 
alized by using "ghost particles" (iTakeda et al. I [1994). 
The time evolution of Rh{t) is obtained by the results 
of the ID simulations. At the four surface areas {y = 
±a;tan(?b, z — ±a;tan^b) in the pyramid-shaped cal- 
culation domain (see Figure [2]) , the rotational periodic 
boundary condition is imposed. 

Since the gravitational force is a long-range force, we 
have to take into account the particles in the whole solid 
angle. However, we only have information in the part 
of the solid angle as shown in Section 12.31 In this pa- 
per, the gravitational force is evaluated in the following 
method. Figure |3] shows the schematic picture of the 
gas sphere viewed from the a;-direction for the case with 
6'b = 27r/10. The latitude and longitude lines are plotted 
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at interval of 0h- The region is surrounded by the thick 
sohd hues represents the calculation domain (see Figure 
[2]). We rotate the calculation domain i6h degrees in the 
^^-direction, where i = —■k/{29\,),...,tt/{29\^). At each 
rotation step in the ^-direction, we rotate the domain 
fc^b degrees in the 0-direction, where k = 1, . . . , 27r/0b- 
After that, the information in the whole solid angle can 
be obtained. However, for i ^ 0, the rotated domains 
in the ^-direction have overlaps. This comes from the 
fact that the pyramid cannot fill 3D space. For example, 
the dashed regions in Figure |3] correspond to the rotated 
domains in the ^-direction. We remove the overlapped 
SPH particles with | tan~^(y/a;)| > 0b from the rotated 
domains. In the simulation, we rotate not SPH particles 
but the tree structure, and calculate the gravitational 
force from the rotated tree structures at each rotation 
step. 

By this method, supposed particle distribution in the 
whole solid angular 47r is not cyclic completely. However, 
the effect of the approximation is negligible since Oh is 
very small in this paper. For i = ±1 (near the equatorial 
plane) , the number of removed SPH particle is negligible 
compared with the total number. On the other hand, for 
large \i\ (near the poles), the fraction of removed SPH 
particles becomes large. However, since the rotated do- 
main is very far from the calculation domain, the detail 
particle distribution is not important. 

In this paper, we neglect the gravitational force outside 
the SF to prevent the ambient gas from collapsing toward 
the center. We find the radius i?grav of the most distant 
particle from the center th at has density p > p th_, where 
pth is a threshold density (jBisbas et al 112003 ). In our 
simulations, we adopt pth = l-lpE- Only SPH particles 
within i?grav are assumed to feel the gravitational force 
that is calculated above method. This means that there 
is the discontinuity in the gravitational force at the i?grav 
In the Appendix, we investigate the effect of the discon- 
tinuity, and confirm that it does not influence our results 
as long as -Rgrav is inside the shock transition layer. The 
main reason is that the pressure suddenly changes in the 
shock transition layer, suggesting that the pressure gra- 
dient is much larger than the jump in the gravitational 
force. 



2.6. Measures of Fluctuations 

In order to evaluate growth of perturbations in the 
SPH calculations, we introduce indicators for fluctua- 
tions of density and the position of the CD. The solid 
angle fib is divided into cells, and the direction of the 
i-th cell is deflned by the unit vector {1 < i < N^). 
The radial proflles of the physical quantities along 
are obtained from the SPH calculations. Here, we deflne 
the angle-dependent maximum density Pmax(ni) along 
rii. The position of the CD along is determined as the 
position where the temperature coincides with a criti- 
cal value that is set to \f2Tc- The speciflc choice of the 
critical value is not important in our results. We aver- 
age pmax(ni) and -RcD(ni) over all directions. The mean 
value of Q = (pmax, ^cd) is defined by 




Fig. 3. — Schematic picture viewed from the a;-direction for 
the case with = 27r/10. The latitude and longitude lines are 
plotted at interval of d\^. The region is surrounded by the thick 
solid lines represents the calculation domain (see Figure [2}. The 
dashed regions correspond to rotated domains in the ^-direction. 

As the indicators of the fluctuations, we evaluate the 
dispersions of these quantities 




2.7. Unperturbed Evolution 

As a test calculation, we present the evolution of the 
shell without initial perturbations in the HM model. Fig- 
ure |4] shows the time evolution of (Rcd) (Equation (fT2|) ) 
evaluated from the SPH simulation (the open circles). 
For comparison, we plot the result of the ID simulation 
by the solid line. One can see that the SPH simulation 
can reproduce the results of the ID calculation very well. 

Figures [S] show snapshots of density (the upper panel) 
and pressure (the lower panel) profiles for t/to = 0.5, 
0.7, 1.0, and 1.26. The abscissae indicate the distance 
from the CD at each time. Only SPH particles in 
\y\/Ro < 0.05 and \z\/Ro < 0.05 are plotted by the dots. 
For comparison, the density profiles in the ID simulation 
are superimposed by the thick gray lines in the upper 
panel. It is seen that the result of the SPH calculation 
agrees with that of the ID simulation very well. Our 3D 
simulation clearly represents the structure and evolution 
of the shell although it is very thin. 

Note that the pressure profiles in the lower panel of 
Figure [5] are smooth at the CD (r — {Rqu) — 0) thanks 
to the modified Godunov SPH shown in Section [2. 1.1 1 In 
the standard SPH, a wiggle in pressure profile appears 
at the CD. In order to see whether the pressure profile 
is smooth even when perturbation is added, the density 
and pressure profiles for t — 1.0 and I — 52 are plotted 
in Figure [HI The detailed investigation of the results 
with perturbation is shown in Section |4l The abscissa 
is the distance from (Rcd)- We plot SPH particles in 
(j)o - A(j) < (j) < (j)o + '^4> and |z|/_Ro < 0.05, where 
A0 = O.OOl^b- Figures[S{a) and (b) correspond to (j) = 
where 6p has the maximum value and 0o = 0b/2 where 
6p has the minimum value, respectively. One can clearly 
see that the pressure profile is smooth at the CD even 
with perturbation. 
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-Rst(pc) 


-Ro(pc) 


To(Myr) 


Mo(A//q) 


High mass (HM) 


41 


0.56 


5.5 


1.6 


4.1 X 10^ 


Low mass (LM) 


19 


0.25 


3.9 


1.6 


1.46 X 10^ 



" Mass of the star 



TABLE 1 

Parameters of Our Models 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 



t/to 



Fig. 4. — Time evolution of (iJcD ) • The open circles indicate the 
result of the SPH calculation. The solid line indicates the result of 
the ID calculation. 

3. PREDICTION BY LINEAR ANALYSIS 

llwasaki et al. I ()201lD performed a linear analysis tak- 
ing into account the thickness of the shell as in Figure 
[5] and imposing the approximate SF and the CD bound- 
ary conditions. They neglected evolutionary effect, i.e., 
expansion and mass accretion through the SF. At each 
instant of time, they solved eigen-value problem and ob- 
tained the growth rate r{l,t) as a function of the time 
and the angular wavenumber, where F is the imaginary 
part of t) in Paper I. The angular wavenumber / can 
be expressed as 27ri?s/fc by using the wavenumber k in 
Paper I, where Rs is the shell radius. 

3.1. Time Evolution and Scaling Law of Density 
Profile 

The density profile of the decelerating shell is asym- 
metric with respect to the density peak, llwasaki et al. I 
(|2011f ) developed a semi-analytic method for deriving 
the time evolution of the density profile by assuming 
that the shell reaches the hydrostatic equilibrium among 
the pressure gradient, the inertia force owing to the 
deceleratio n, and the radial grayitatio n a.1 force toward 
the ce nter (jWhitworth fc Francis II2Q02I ). llwasaki et al. I 
(|2011l ) found that the time evolution of the density 
profile can be divided into three evolutionary phases, 
deceleration-dominated, intermediate, and self-gravity- 
dominated phases (see Figure 2 in Paper I) . In the early 
phase (deceleration-dominated phase), the density peak 
is at the SF because the inertia force is larger than the 
gravitational force. As the shell expands, the inertia force 
decreases while the gravitational force increases. Thus, 
the density peak moves from the SF toward the CD as 
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Fig. 5. — Snapshots of density (the upper panel) and pressure 
(the lower panel) profiles for t/to = 0.5, 0.7, 1.0, and 1.26. Only 
SPH particles in \y\/Ro < 0.05 and \z\/Ro < 0.05 are plotted 
by the dots. The abscissae are the distance from the CE). The 
thick gray lines in the upper panel indicate the results of the ID 
simulation. 

shown in Figure |SJ In the intermediate phase, the den- 
sity peak is closer to the SF than the CD. In the later 
phase (self-gravity-dominated phase), the density peak is 
closer to the CD than the SF. 

It is also found that the density profiles for various sets 
of (tie, Quv) E^re characterized by a single parameter, 
that is the typical Mach number, 

4 i?o _ „l/7 „-l/2 -1/14 , 
-'^O - 3— - 7 guv.49 -'c,10 "e,3 ' (14) 

where T^^w = TjlO K. 

They found the development of the GI is strongly in- 
fluenced by two effects, asymmetry of the density profile 
and boundary conditions. 

3.2. Asymmetry of the Density Profile 

The asymmetry of the density profile strongly influ- 
ences the GI, especially, in the self-gravity-dominated 
phase (for example, the shell at t/to = 1-26 in FigurejS]). 
In this later phase, the distance from the density peak 
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Fig. 6. — Snapshots of the density and pressure profiles for t/f = 
1.0 and I = 52. The abscissae are the distance from the average 
position of the CD. SPH particles in 0o — ^0 < < 0o + ^<f> a-nd 
\z\/Ro < 0.05 are plotted, where (po =(a)0, {h)0^,/2, and A<j> = 
0.0010b. 

to the SF is larger than Hq, where Hq — Cg/ \/2ttGpoo is 
the scale height of the shell and poo is the peak den- 
sity. On the other hand, the distance from the den- 
sity peak to the CD is smaller than Hq. The gas tends 
to collect toward the density peak because the unper- 
turbed gravitational potential has the minimum value 
there. The gas near the SF can collapse toward the den- 
sity peak because the sound wave cannot travel from 
the density peak to the SF. This mode is "compress- 
ible mode". On the other hand, the gas near the CD 
cannot collapse to the peak. However, the GI can pro- 
ceed even there through the deformation of the CD that 
makes the gravitational potential deeper. This mode is 
"incompressible mode" (jElmegreen fc Elmegreeni 119781 : 
iLubow k Pringle I [1991 . Therefore, the GI of the sheh 
in the later phase has properties of compressible and in- 
compressible inodeS; 

Moreover. Ilwasaki et al. I ()2011t ) found that the growth 
rate is enhanced compared with symmetric case through 
the combination of the GI and the Rayleigh- Taylor in- 
stability. 

3.3. Boundary Conditions 

Expanding shell is confined by the SF on the lead- 
ing surface and by the CD on th e traiU ng surface. 
Dispersion relations of lElmegreen I (| 19941 . E94) and 
iWhitworth et al."l (|1994[ ) are essentially based on the dis- 
persion relation of the layer confined by the SF on both 



surfaces. The shock-confined layer is stabilized by the 
tangential flow generated by obliqueness of the SF com- 
pared with the layer confined by thermal pressure of hot 
gas. Therefore, the linear analysis in Paper I taking into 
account SF+CD boundary conditions gives larger growth 
rate compared with E94. 

In summary, as described in Sections l3.2l and [5751 Paper 
I predicts that the GI begins to grow earlier than the 
prediction from E94. The development of the GI in the 
later phase is accompanied by the significant deformation 
of the CD. 

4. RESULTS 
4.1. The Case with High-mass Central Star 

First, we present the results of the HM model with 
perturbations. We calculate four runs for the case with 
I = 26, 52, 78, and 156. Figures [H; a), (b), (c), and (d) 
represent the time evolution of the perturbation ampli- 
tude for the wavenumber I = 26, 52, 78, and 156, respec- 
tively. The thick solid and the thick dashed lines corre- 
spond to the density perturbation Apmax/ (Pmax) and the 
deformation of the CD 30 x Ai?cD , respectively. In each 
panel, Apmax/ (pmax) and Ai?cD grow in the similar way. 
Perturbations with / = 52 — 78 grow with the largest 
growth rate. For larger angular wavenumber I = 156 
and smaller angular wavenumber I = 26, they grow more 
slowly. The prediction from E94 is also superimposed by 
the thin dashed line in each panel. It is found that in 
our 3D results, perturbations begin to grow earlier than 
the prediction of E94. 

Let us compare the results of SPH simulations with the 
linear analysis in Paper I. However, Paper I neglected the 
evolutionary effects, i.e., expansion of the shell and the 
mass accretion through the SF. To include these effects 
approximately, we modify the instantaneous growth rate 
from r to Fevo as follows: 

where is the velocity of the shell. This modification is 
inspired by E94. The time evolution of Rs{t) and Vs{t) is 
determined by the thin-shell model shown in Section 2 in 
Paper I. One can see that the terms of Vg/Rs correspond 
to evolutionary effect of the shell that stabilizes the GI. 
The prediction based on Paper I is superimposed by the 
thin solid line in each panel of Figure [T) It is evaluated 
by exp [J^ Tcm{l,t')dt'] . From Figure [3 one can see that 
the linear analysis of Paper I describes the growth of 
perturbations very well. 

Figure [5] shows the time sequence of cross sections 
through the 2 = plane for I — 52. The color scale indi- 
cates density, and the arrows represent relative velocity 
to the fluid with the maximum density. The abscissa 
and the ordinate axes correspond to the azimuthal an- 
gle (j), and the distance from the CD divided by (Rcd), 
respectively. Figure [Sfa) represents the initial condition 
where the shell around cj) = is displaced to the neg- 
ative radial direction. Figure ^h) shows that the tan- 
gential flow is generated by the obliqueness of the SF, 
and collects the gas around = 0. As a result, en- 
hanced pressure pushes the SF outward in the radial 
direction (see Figure El^c)). The gravitational potential 
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around = becomes deep, leading to gas accumula- 
tion (see Figure [HKd)). This mode of the GI is similar to 
the incompressible mode in th e sense that the deforma- 
tion develops at the boundary (jElmegreen fc Elmegreen I 
[Tots'; 'Lubow & PringleJll993|). In the later phase (Fig- 
ure HJc)), the CD highly deforms while the SF is almost 
flat. This is consistent with prediction of Paper I (see 
Section [S]). Figure |8l[e) is very similar to Figure 11 in 
Paper I that shows the eigenfunction in parameters of 
the HM model at t/to = 1.3 for I = 52. 

Figure [9] is the same as Figure |8] but for larger 
wavenumber case, I = 156. The initial state is shown 
in Figure inija). Because of large wavenumber, the tan- 
gential flow collects the gas more quickly in Figure [9Kb) 
that is similar to Figure [UJc). In this moment, since the 
SF at = deforms to the positive radial direction, 
the gas moves away from (/> = by the tangential flow. 
As a result, at t/to = 0.8 (Figure IlKc)), the SF and the 
CD are almost flat while the velocity field exists inside 
the shell. The density distribution at t/^o = 1-0 (Fig- 
ure [Qjd)) has the opposite phase compared to the initial 
state. In this moment, the perturbation begins to grow 
but the growth rate is smaller than that for 1 — ^2. The 
CD deforms largely at the final state in both of Figures 
[5] and ini In this sense, the behavior of the GI is similar 
to that for / — 52. 
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Fig. 7. — Time evolution of perturbations for the HM model. 
Each panel corresponds to Z = (a) 26, (b) 52, (c) 78, and (d) 
156. The thick solid and the thick dashed lines represent the den- 
sity perturbation Apniax/(pmax) and the deformation of the CD 
30 X Ai?cDi respectively. The thin solid lines correspond to the 
prediction from Pa per I. The th in dashed lines correspond to the 
prediction based on lElmegreen I fl994) . 



4.2. The Case with Low Mass Central Star 

Next, we present results of the LM model. Because the 
central star is less energetic, the Mach number of the ex- 
panding shell is smaller than that in the HM model, sug- 
gesting that the shell is geometrically thicker. Therefore, 
perturbations are expected to grow more slowly than the 
HM model. This feature is confirmed in Figure [TU] that is 
similar to Figure [7| The most unstable mode is found in 
smaller wavenumber I 28 and the growth rate is smaller 



Fig. 8. — Cross sections through the z = plane for I = 52. 
Each panel shows the shell at t/to = (a)0.4, (b)0.6, (c)0.8, (d)l.O, 
and (d)1.3, respectively. The color scale represents the density, and 
arrows are velocities relative to the fluid at the density peak. The 
abscissa and the ordinate axes correspond to the azimuthal angle 
(j> = tan~^ {y /x), and i?/(i?cD) ~ li respectively. 



than the HM model. We found that perturbations be- 
gin to grow earlier than the prediction of E94 also in this 
case. The linear analysis in Paper I describes the GI very 
well. Figure [TT] shows the cross section through z = for 
^ = 56 at t/to = 1.5. The significant deformation of the 
CD is seen as in the HM model. 

5. DISCUSSION 

5.1. Crowth Rate of Gravitational Instability 

We presented only two examples of the HM and the 
LM models in Section In this section, we predict the 
GI in the shells with various parameters (tie, QuVi 7c) 
by using the results of the linear analysis presented in 
Paper I. As shown in Figures [7] and [TUl it is found that 
the results of the SPH simulations are well described by 
the growth rate, Fcvo- Using this growth rate, one can 
predict the time tbgn when the GI begins to grow by the 
condition rovo(^bgn, ^bgn) — 0, where ^bgn is the angular 
wavenumber of the most unstable mode at t = tbgn- This 
condition is rewritten as 



V877" — r(^bgn, ibgn), 



(16) 



where the left-hand side corresponds to the evolution- 
ary rate owing to expansion and mass accretion, the 
right-hand side corresponds to the growth rate of the GI. 
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Fig. 9. — Same as Figure[8]but for / = 156. 




Fig. 10. — The same as Figure [7] but for the LM model. Each 
panel corresponds to I =(a)14, (b)28, {c)56, and {d)84. 



At early phase, the evolutionary rate is larger than the 
growth rate of the GI, indicating that the shell is stable. 
As the shell expands, Vs/Rs decreases while F increases. 
Therefore, the GI begins to grow when the growth rate 
of the GI is larger than the evolutionary rate. The radius 
of the shell at t = tbgn is defined by i?s.bgn- 
Figure fT2l a). (b), (c), and (d) show contour maps of 




Fig. 11. — Cross section through 2 = plane for Z = 56 at 
t/to = 1.5 for the LM model. The color scale represents the density, 
and the arrows show velocities relative to the fluid at the density 
peak. The abscissa and ordinate axes are the same as Figure |8] 



bgn 



t 



|7r(A- 



bgn 



A, 
/2)2 



27ri?s_bgn/^bgn, and Mbgn = S(t = 
in the (rig, Quv) plane derived by us- 



ing the linear analysis, respectively. Here, Tc is as- 
sumed to be 10 K. From Figures [121 one can see that 
the ibgn, ^s.bgn, Abgn, and Mbgn decrease with increas- 
ing ng. This can be understood by scaling of the typical 
number density of the shell, nsh = ueMq that corre- 
sponds to the shocked gas density. From Equation (fT4| . 

we have Ugh oc n'^'^ ■ Therefore, for larger ue, the shell 
becomes denser and unstable earlier (Figure [T^a)). and 
the most unstable scale is smaller fFiguresfT2lc) and (d)). 
Next, let us see the dependence of ibgn, -Rs,bgn, Abgn, 
and Mbgn on Quv- Since the central star is more ener- 
getic for larger Quv, the Mach number of the expanding 

shell (m oc Quv) larger, indicating that the shell 

is denser and is more unstable. Therefore, ibgn, Abgn, 
and Mbgn decreases with increasing Quv- On the other 
hand, i?s,bgn increases with Quv because the expansion 
velocity increases with Quv- 

From Figures I12| since the contour lines are almost 
straight in the plane (logj^Q tie, log^Q Quv), quantities are 
expected to have power-law dependence on Quv and tie, 
or equivalently on A4q. In Figure we plot tbgn/^o, 
Rs.hgn/Ro, and Zbgn as a function of Mq for various pa- 
rameters that correspond to the open circles in Figurc[T2l 
The open circles correspond to the results of the linear 
analysis. From Figures [THl one can see that the results of 
the linear analysis are well described by the analytic for- 
mulae, tbg„/to = 1.62^^o°■^^ Rs^jRo = 1.2Xo"■3^^ 
and Zbgn = 3.2A^Q ^ that are plotted by the solid lines. 
Using Equation (|14p , the analytic formulae are rewritten 
as 



tbgn = 0.6 Q 



-0.11 rpO.375 



UV,49 -'cao 



-0.45 

*E,3 



Myr, 



Rs 



O.'i (^UV,49 -'c.lO "E,3 



0.19 „-0.54 



and 



^bgn 



54 Q 



0.21 



UV,49 -^c,10 



-0.75 „-0.11 



(17) 
(18) 



(19) 



The wavelength Abgn and the mass A/bgn that correspond 
to Zbgn are given by 

Abgn = 0.4 Q^"^%' T,"f 4 „-o.44 (20) 



and 



Mb 



o t- ^-0.16 ^2.06 „-0.42 T>/r 
^■O '^UV,49 ^c,10 '*E,3 ^"0' 



(21) 
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Fig. 12. — Contour maps of (a)ii3gn, (bVRg^bgm (c)Abgni and (d)Mbgn in the (logj^gnE, logj^Q Quv) plane. The open circles correspond 
to the parameters that are used in Figure fT3l 




Fig. 13. — Dependence of (a)tbgii/to, (b)Hs,bgn/^Oi and (c)iniax 
on the typical Mach number Mo- The open circles correspond to 
the results of linear analysis. The solid lines indicate the analytic 
formulae. 



iWhitworth et al.~l ()1994 W94) also derived similar for- 
mulae by using the thin-shell linear analysis similar to 
E94. The dependence of Equations (IT71) - (PT]) on param- 
eters QuM, Tc, and shows close agreement with that 
in W94 although the numerical factors are different. Our 
results show that the GI begins to grow earlier and as a 
result, the fragment mass is smaller than their results. 
In detail, ibgn, -Rbgn, and Abgn are roughly half of those 
in W94, and Mbg„ is roughly 1/8 of that in W94. This is 
because the stabilization effect of the SF does not work 
on the trailing surface that is the CD as shown in Paper 
I. As pointed out in W94, the properties of fragments 
are insensitive to QuVi and mainly depend on tie and 
Tc- Indeed, in Equations and (PT|) . the dependence 
of Abgn and Mbgn on tie are close to the Jeans length and 

the Jeans mass of the shell (Aj, Mj) cx n^^^'^ oc n-^^^'^ . 

5.2. Comparison with Previous Studies 

lElmegreen I ([1989) investigated the GI in a decelerating 
shocked layer. He numerically integrated perturbation 
equations with respect to time by taking into account 
time evolution of the layer and by averaging physical 
quantities over the thickness. He has taken into account 
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SF and CD boundary conditions on the SF and CD, 
respectively. He called the boundary effect on the CD 
"pinching force" . In his linear analysis, the destabiliz- 
ing effect of the pinching force does not work efficiently. 
One reason is that he used Equation (3) in Paper I as 
the pressure of the HII region although the geometry 
is plane-parallel. Therefore, the pressure at the CD is 
underestimated, leading to underestimate the pinching 
force of the CD. The other reason is the geometry of 
the gravitational potential. Asymmetry of density pro- 
file of spherical shell is qualitatively different from that 
of plane-parallel layer. In the plane-parallel geometry, 
the static symmetric layer peaks at the mid-plane. If the 
layer decelerates, the density peak is always closer to the 
SF than the CD. Therefore, the tangential ffow behind 
SF erodes the density perturbation more efficiently. On 
the other hand, for the case with shells, the peaks can 
be closer to the CD than the SF (see Section [Q]) . From 
these reasons, he derived lower growth rates than our 
r esults. 

iDale et all ()2009f ) simulated the gravitational fragmen- 
tation of expanding shells confined by temporary con- 
stant thermal pressure of hot rarefied gases on both sides. 
In their calculation, the motion of the shell is determined 
so that the pressure at the leading surface is the same 
as that at the trailing surface. Therefore, the density 
peak is always around the mid-plane of the shell, and 
the density profile is almost symmetric. In their calcu- 
lation, the column density decreases with time because 
the shell expands keeping the mass fixed. Therefore, the 
pressures at the boundaries approach to the peak pres- 
sure. They found that the confining pressure accelerates 
fragmentation in the later phase, and describe this effect 
as "pressure-assisted" gravitational fragmentation. This 
mode is the same as the incompressible mode. Recently, 
taking into acc ount the effect of the t hickness of the shell 
approximatelv. IWiinsch et al. I (|2010[ ) established a semi- 
analytic linear analysis that explains results of their SPH 
simulations. Different from them, in Paper I, we have 
taken into account the effect of SF-I-CD boundary con- 
dition approximately. 

5.3. Other Instabilities 

We discuss other potentially important effects that are 
not analyzed in this paper. 

5.3.1. Vishniac Instability 

iVishniad ()1983D found a hydrodynamical overstability 
on decelerating shells confined by the SF on the leading 
surface and by the CD on the trailing surface by using 
thin-shell linear analysis. However, the Vishniac insta- 
bility (VI) did not appear to occur in our simulations. 
This can be explained by the stabilizing effect by expan- 
sion. iVishniac fc Rvu I ()1989[ ) obtained an analytic dis- 
persion relation of the VI of the expanding shell without 
self-gravity by taking into account the thickness approx- 
imately. In Section 6 in Paper I, they discussed the VI 
by applying the analytic dispersion relation to the shell 
driven by the HII region. They found that the VI is ex- 
pected not to be important in the later phase {t/to > 0.5) 
since its growth rate is comparable to the expansion rate 
of the shell radius. Therefore, the VI was not found in 
our simulation. However, they also found that before the 



beginning time of our simulations {t/to < 0.4), the VI is 
expected to grow rapidly since the Mach number is large 
(see Figure 13 in Paper I). The perturbations quickly 
reach the nonlin ear regime and saturate as a subsonic 
transverse flow (|Mac Low fc Norman I [19931 ). As a re- 
sults, a small scale subsonic turbulence whose angular 
wavelength of I — 10^ ^ 10^ is expected to be induced. 
The consequence of VI may correspond to the increase 
of the initial perturbations in our 3D simulation. 

5.3.2. Thermal Instability 

In this paper, we assume that the fluid evolves keep- 
ing its temperature fixed. However, in real interstellar 
medium (ISM) , shock-heated gas cools via radiative cool- 
ing. It is well known that the cooling ISM is often ther- 
mally unstable. In such a case, during the cooling con- 
densation, the layer is expec ted to fragment into dens e 
clouds with various scales f Iwasaki fc Tsuribe 1 12009^. 
Koyama fc Inutsuka (2002) investigated the propaga- 
tion of a shock wave into the warm neutral medium. 
They found that cold cloudlets move with supersonic ve- 
locity dispersion in surrounding warm gas in the post- 
shock region. The velocity dispersion is larger than the 
sound speed of cloudlets while it is smaller than the 
sound speed of surrounding warm gas. Since the shocked 
molecular cloud is also thermally unstable, the thermal 
instability drives supersonic turbulence inside the shell, 
and cold molecular cloudlets move with the supersonic 
velocity dispersion. This situation is quite different from 
the isothermal gas that has been considered in this pa- 
per. The supersonic turbulence is expected to be im- 
portant in contrast to the subsonic turbulence driven 
by the VI, and considerably modify the evolution. The 
mass of the cold molecular cloudlets increases with time 
through the accretion of surrounding gas and the coales- 
cence between cloudlets. The supersonic turbulence may 
provide effective turbulent pressure, while coalescence of 
cold cloudlets may decrease turbulent pressure. We will 
address the effect of the thermal instability on the GI of 
the expanding shells in our next work. 

5.4. Ionization Front 

In the paper, we did not solve the radiative transfer 
equation of the ionizing photons. Instead, we assumed 
that the trailing surface is the CD instead of the ioniza- 
tion front (IF). There are two effects of the ionization 
on the Gl. One is that a fraction of the shell becomes 
ionized gas as HII region expands. Therefore, our simu- 
lations overestimate the column density of the shell. The 
swept-up mass by the SF is given by 

Msw = yPEi^s- (22) 
On the other hand, the ionized mass is given by 

Mion = YPii?^ (23) 

where pi is the density of the ionized gas, and we neglect 
the difference between radii of the IF and SF. The real 
shell mass is given by Mgw — Mion- Inside the HII region, 
the balance between the ionization and the recombina- 
tion is approximately established. Therefore, the density 
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of the HII region becomes 



uv 



1/2 



PE 



R. 



•ST 



3/2 



(24) 



where m is the mean weight of the gas particles and we 
use Equation ([T]). Therefore, from Equations (|22)) - ((24)) . 
the ratio of A/ion to Afgw is 



7^ 



Ms, 



(25) 



If 7?. <C 1, the effect of the ionization is negligible. In Fig- 
ure 1 in Paper I (also see Figure^), the Stromgren radius 
is as large as 0.1i?o in various parameters (Quv, "-e)- 
At the beginning time of our simulations, t — OAIq 
{Rs ~ 0.6i?o), 7^ - (0.1/0.6)3/2 = 0.07. As the shell ex- 
pands, n decreases. At t = to, 7^ - (O.l/lO)^/^ ~ 0.03. 
Therefore, the contribution of the ionization is limited to 
only several percent. 

The other is that the IF induces an instability on 
the shell. The D-type IF is analogous to a combus- 
tion front that is unstable to corrugational deform ations 
(|Landau fc Lifshitz I [l987f ). iVandervoortl ([1962') found 
that the D-type IF is unstable. In this instability, the 
most unstable scale is comparable to the t hickne ss of the 
IF. Moreover, iGarcia-Segura fc Franco I ()1995[ ) showed 
that the VI is strongly modified by the IF. The shell 
rapidly fragments and finger-like structures are gener- 
ated in their two-dimensional simulations. Note, how- 
ever, that they assumed the power-law density distribu- 
tion oc r"™ of the ambient gas with the uniform density 
core in the center with 0.2 pc, where parameter w is set 
to 2. The expansion of the IF depends on the power- 
law index w in the ID model, li w > 1.5, even if the 
SF is emerged in front of the IF in the uniform density 
core, the IF quickly gets ahead of the SF, and eventually 
the whole of the shell and t he ambient gas are ionized. 
Moreover, iHosokawal (|2007[ ) investigated the expansion 
of the IF and the dissociation front, and showed that the 
star formation is expected to be suppressed when w > 1. 
This is because the column density decreases as the shell 
expands and the FUV photon easily escapes in front of 
the shock. Therefore, model bv lGarcia-Segura fc Francol 
([1995) corresponds to the case where the triggered star 
formation is not simply expected. We need to investi- 
gate the expansion of the IF for the case with the shal- 
lower density profile w < 1 in detail using the radia- 
tive multi-dimensional sim ulations taking in to account 
the self-gravity. Moreover, iWilliamsl (|1999f ) discovered 
shadowing instability in the R-type IF before emerging 
of the SF. This instability may also disturb the shell in 
the very early phase of evolution. 

5.5. Magnetic Field 

It is well known that the magnetic field is important 
in the dyna mics of the ISM alt hough it is neglected in 
this paper. iNagai et al. I (|1998?) investigated the GI of 
a pressure-confined isothermal layer threaded by uni- 
form magnetic field by using linear analysis. They found 
that the magnetic field cannot stabilize the GI. More- 
over, they found the layer fragment in filamentary clouds 



whose longitudinal axis is parallel or perpendicular to the 
magnetic field depending on the thickness of the layer. 
However, their analysis is restricted to the static mag- 
netized layer. It is important to investigate stability of 
magnetized shocked layer or shell in dynamical modeling. 

6. SUMMARY 

In this paper, we have investigated the GI of the ex- 
panding shell using 3D numerical simulation with high 
resolution. We summarize our results as follows: 

1. The GI begins to grow earlier than the prediction 
from the line ar analysis based on the thin-shell ap- 
proxi mation (jElmegreen I Il994t iWhitworth et al."1 
119941 ). During the development of the GI, the con- 
tact discontinuity highly deforms while the shock 
front remains almost flat. These all features are 
consistent with the prediction from llwasaki et alTI 

dmiD. 

2. The GI is expected to begin to grow at an epoch 
when the growth rate of the GI becomes larger than 
the evolutionary rate owing to the expansion and 
the mass accretion (see Equation (IT6l) ) Using the 
results of the linear analysis, we have derived use- 
ful approximate analytic formulae (Equations ()17|) - 
([2T|l ) for the fragment scale and the epoch when the 
GI starts. 

This paper and Paper I revisit the fragmentation of the 
expanding shells by using both the 3D simulation and the 
linear analysis. Since the gas is subject to heating owing 
to far-UV photo n, the gas temperature is larger than 10 
K (|Hosokawa &: Inuts uka 2005, 2006). When Tc = 30 K 
(50 K), the GI with mass scale of Mbgn 34Mq (96 Af©) 
begins to grow at t^gn 0.9 Myr (1.1 Myr) (Quv = 10*^ 
s-i, riE = 10^ cm"3) from Equations and ^l^j. The 
fragment mass is expected to be larger than A/bgn be- 
cause it increases through the gas accretion. Moreover, 
perturbations with larger scale than the most unstable 
mode is also unstable. Therefore, the larger scale pertur- 
bations gradually grows, and small scale fragments may 
coalesce into larger one, and the total number of frag- 
ments may decrease. The predicted masses are roughly 
comparable to (or are slightly smaller than) observational 
masses of dense cores (dozens ^ hundreds M ©) (e.g., 
iDeharveng et al. I l2003HZavagno et al. Il2006f ) although 
our model is based on very idealized situation, such as 
uniform ambient gas and isothermal equation of state. 
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APPENDIX 

EFFECT OF DISCONTINUITY IN GRAVITATIONAL FORCE 

In this paper, the gravitational force is neglected outside i?grav (see Section [275]) . We investigate how the sudden 
change in the gravitational force at i?grav influences our results. To see this, we focus on the time evolution of the 
displacement of the SF ARgp whose definition is the same as Ai?cD (see Section 12. 6|) . When Ai?sF is evaluated, 
the position of the SF along is determined as the density coincides with l.lpE- Figure [TWa) shows that the time 
evolution of AiJgp for / — 52 in the HM model. The solid line and the open circles correspond to pth = 1-lpE and 
2pE, respectively. One can see that they agree very well. Therefore, our results are insensitive to the position of the 
discontinuity of the gravitational force. This is because the pressure gradient is much larger than the gravitational force 
at the shock transition layer. This is clearly seen in Figure [TU^b) that shows the pressure gradient, the gravitational 
force and the density at t — tg for / = 52 in the simulation unit. The SPH particles around cj) — is plotted. The 
threshold density is set to l.lpE- One can see that the gravitational force has the jump around (r— {Rc))/Ro ^ 0.013 
that corresponds to i?grav in Figure [TW b). From Figure [T^ b). the discontinuity in the gravitational force is much 
smaller than the pressure gradient because the pressure suddenly changes at the shock transition layer. 

Therefore, it is concluded that the existence of the discontinuity in the gravitational force is not important as long 
as the discontinuity is inside the shock transition layer. 




Fig. 14. — a) Time evolution of the displacement of the SF. The solid line and the open circles correspond to = l.lpE a-nd 2pEi 
respectively, b) Snapshots of the pressure gradient (the light gray circles), and the gravitational force (the dark gray circles) along the 
r-direction in the simulation unit. The density distribution is superimposed by the filled circles. The abscissae are the distance from the 
CD. 



